A logarithmic contribution to the density of states of rectangular Andreev billiards 
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We demonstrate that the exact quantum mechanical calculations are in good agreement with 
the semiclassical predictions for rectangular Andreev billiards and therefore for a large number of 
open channels it is sufficient to investigate the Bohr-Sommerfeld approximation of the density of 
states. We present exact calculations of the classical path length distribution P(s) which is a non- 
differentiable function of s, but whose integral is a smooth function with logarithmically dependent 
asymptotic behavior. Consequently, the density of states of rectangular Andreev billiards has two 
contributions on the scale of the Thouless energy: one which is well-known and it is proportional 
to the energy, and the other which shows a logarithmic energy dependence. It is shown that the 
prefactors of both contributions depend on the geometry of the billiards but they have universal 
limiting values when the width of the superconductor tends to zero. 
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When a normal-metaUic dot is placed in contact with a 
superconductor, the low-energy density of states n[E) of 
the resulting Andreev billiard is strongly modified com- 
pared with the normal state 0. The energy depen- 
dence of n{E) is highly non-trivial and provides a testing 
ground for current understanding of proximity effects of 
hybrid superconducting nanostructures and for theoreti- 
cal tools such as semiclassical theory. Recently a number 
of conflicting results have been obtained for the density 
of states (DOS) of ballistic Andreev billiards (AB). One 
example is shown in figure 1, which in the separable limit 
W = a was studied long ago in 0, where it was found 
that for small E, n{E) is proportional to E. In contrast 
for W < a, the problem has only recently been studied 
B II II IE 13 • In Ref. it was found that n{E) = vE, 
where v is a universal constant, independent of E and 
W. In contrast 0] predicted different values for v. 

In this paper we resolve this discrepancy by showing 
that neither of these results is strictly correct in the limit 
E ^ and instead predict that n{E) diverges logarithmi- 
cally. This generic logarithmic contribution is significant 
for energies less than or of order the Thouless energy. 

To address this problem, we use both exact and semi- 
classical techniques to analyze the AB of figure 1. The 
exact calculation starts with a ballistic two dimensional 
normal dot of area A, described by a scattering matrix 
So [E) and with a mean level spacing for the isolated nor- 
mal system 5 = at the Fermi energy Ep. Then, a 
superconductor of width W and bulk order parameter A 
is placed in contact with such a billiard. The number 
of open channels in the S region is the integer part of 
M = ''^J^ , and the energy levels of AB are the posi- 
tive eigenvalues E (measured from the Fermi energy) of 
the Bogoliubov-de Gennes equation^. A secular equa- 
tion of the exact energy levels of AB, in terms of the 



scattering matrix Sq{E) of the normal region, is derived 
by matching the wave functions at the interface of the 
normal-superconductor (N-S) systems. The energy lev- 
els are exact in the sense that no Andreev approximation 
{A/Ep <C 1 and quasi-particles whose incident/reflected 
directions are approximately perpendicular to the N-S 
interface) ^ is assumed. 

We also give the semiclassical Bohr-Sommerfeld ap- 
proximation of the density of states (DOS) n{E) ex- 
pressed by the classical return probability P{s) of the 
electron. Our exact and semiclassical calculations are 
applied to the integrable AB shown in Fig.^ We demon- 
strate that the integrated DOS agrees very well with the 
exact calculations. We find that the integral of the re- 
turn probability P{s) has a contribution depending on 
s logarithmically in the asymptotic limit, s —f oo, which 
has to-date been overlooked. As a consequence, the small 
energy dependence of the DOS n{E) has also a logarith- 
mic factor in addition to a contribution which depends 
linearly on the energy E as predicted in Ref. 0, 0. 
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FIG. 1: A normal dot of rectangular shape in contact with a 
superconductor placed at the middle of the edge of the normal 
region. 

In our earlier work we found that the exact energy 
levels of AB with arbitrary shape of normal region are 
the solutions of the following secular equation: 

det(^Im{7Do(E)Dh(E)}^ = 0, where (la) 
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De{E) = Q{E)+K{E)G{E), (lb) 
Dh{E) = [Q{-E)-K{-E)G*{-E)]-\ (Ic) 
G{E) = [l-So{E)][l + So{E)]-\ (Id) 

Here 7 = e~*"'^'^°^(^/'^), Q and K are diagonal matri- 
ces with elements Qnm (E) ~ (5„m In {E) and if ,„„ (i?) = 



kn[E), where g„(-B) = kpx 1 + 



the transverse wavenumbe rs of the ele ctron in the S 

are the trans- 



region and kn{E) = kp^l + ^ ~ jj2 
verse wavenumbers of the electron in the S region when 
A = 0. It is assumed that the Fermi wavenumber, 
kp = \J1mEp jT? is the same in the S and N regions. 
All the matrices are M by M dimensional. All the in- 
formation on the normal region are incorporated in the 
matrix G via the scattering matrix Sq {E) . Note that the 
secular equation Hla() is an extension of that derived for 
box and disk geometries in Ref. 0. 

The density of states in the semiclassical Bohr- 
Sommerfeld approximation is written as Q 



i{E) = M 



/ dsP{s) 
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This expression reduces to that of by Melsen et al. 
Schomerus and BeenakkerQ, Lodder and NazarovQ, 
and Ihra et al. Q, in the limit <C A and when the 
coherence length in the superconductor, ^0 = hvp/A <^ 
d. In that case, the energy dependent phase shift 
— arccos(i?/A), due to Andreev reflections, was approx- 
imated by tt/2. However, our expression is valid with- 
out such an approximation, and essentially improves the 
agreement between the exact and the Bohr-Sommerfeld 
approximation of the DOS (see also Ref. 1^). 

From Eq. 10) one can find a simple expression for the 
integrated DOS N{E) = n{E')dE': 



N{E) = Mj2{'^'F[sn{E)]}, where (3a) 

71=0 

(nTT -f arccos -§) ^ ^ 

sn (E) = ^ ^J^^ Co , (3b) 

F(s) = / P{s')ds', (3c) 



where F{s) is the integrated distribution function of the 
return probability of the electron. Note that P{s) is nor- 
malized to one, i.e., F{oo) = 1. 

To obtain the exact energy levels of the AB shown 
in Fig. ^ we need to calculate the scattering matrix 
So{E) of the N region. Following the same approach 



as Mortensen and co-workers [loj. we obtain 

Sa{E) = qR{X,Q)q^ - I, where 
R{X,Q) = 2(l + g)"'{l 

-2X2 [1 j^x^ + Q(l - X"^)] "^|. 



(4) 



are Here X is a diagonal matrix with elements X„i„ = 
5mn exp{iknd), Q = Q, and g is a [fcpI^/Tr] by [kpa/ir] 



dimension matrix with elements given by the overlap in- 
tegrals defined in Ref. ^3 ([■] stands for the integer part). 
Notice that if = a then g = 1 and the scattering ma- 
trix So{E) = —X^ is a diagonal matrix. In this case, we 
obtain the same secular equation for AB as that in Ref. 
for box geometries. 

In Fig. 121 the exactly (numerically) computed inte- 
grated density of states obtained from Eq. and its 
evaluation in Bohr-Sommerfeld (BS) approximation us- 
ing Eq. are shown for different widths W of the lead 
(for the calculation of P{s) see below). This shows that 
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FIG. 2: The integrated density of states, N{E) from exact 
quantum mechanical calculations (solid lines) and from the 
Bohr-Sommerfeld approximation given by Eq. ^ (dashed 
lines) as functions of E (in units of A) for W = 0.8a (a), 
and W — 0.5a (b). Insets show the enlarged portions ot N{E) 
from exact calculations (solid line), its BS approximation from 
Eq. Q (dotted line) and its asymptotic form from Eq. 18all 
(dashed line). In both cases d = a and the parameters are 
M = 55.5, A/Ef = 0.015. 

the exact calculations and the semiclassical predictions 
are in very good agreement. 

Calculations for P(s) and F{s) start with unfolding the 
trajectory of the electron, i.e. using the fact that the free 
motion of the particle in the billiard is equivalent to its 
flight in a lattice of vertical intervals with length W and 
with lattice constant 2d (a) in the horizontal (vertical) 
direction. In general, one can calculate P{s) by using a 
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large number of trajectories and determining the distri- 
bution of their path lengths. An analytic form of P{s) 
and F{s) can be derived for W > a/2. As an example, 
these functions are plotted for d = a and W ~ a/2 in 
Fig.El Clearly, for arbitrary W, P{s) = for s < 2d and 
for larger s it is a non-differentiable function possessing 
a singularity at s = 2c? and peaks at multiples of 2d. Us- 
ing Eq. one can calculate the DOS and pronounced 
peaks (indeed singularities) arise due to the singularity 
of P{s) at s = 2d. The positions of these singularities 
in the DOS are in perfect agreement with that obtained 
from the general expression derived in Ref. 0. 

The integrated path length distribution F{s) is a 
smooth function of s allowing one to obtain its asymp- 
totic form Fa{s) for s — > oo (see the inset of the top 
panel of Fig.|2Jl. From Eq. Q one can see that the large 
s behavior of P{s) is related to the low energy depen- 
dence of the DOS. Due to the rapid variation of P{s), it 
is difficult to define its asymptote. To avoid this prob- 
lem we first calculate the asymptotic behavior of F{s) 
and the asymptotic form of P{s) is then calculated from 
Pa{s) = dFa{s)/ds (see also the bottom panel of Fig.|3Jl. 
For s — > oo analytical form of the asymptote of F{s) is 
found by considering how particles can travel a large dis- 
tance s without hitting the vertical intervals of length W 
in the unfolded space. Such trajectories lie only in cer- 
tain directions forming corridors with slopes u that are 
multiples of a/{2kd), where k = [W/a]. The probabil- 
ity that a particle runs in such a corridor is proportional 
to 1/(ms^) for large u but |m| < Umax — s/x^-iin, where 
2^min IS a minimal displacement in the horizontal direc- 
tion. Summation of the probabilities for these corridors 
yields a term proportional to In u,„ax- Directions with 
slopes larger than Umax should be treated separately and 
the result contains a term s^^lnccmin- The above terms 
from the two regions sum up to s^^lns, and finally we 
find 



F{s) « Fa{s) = 1 - Cl 



for large s. Here Ci and C2 are functions of d/a and 
W/a. For the case W > a/2 the detailed calculation 
gives analytic expressions also for Ci and C2: 
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(6b) 



(6c) 



Here 7 is the Euler-Mascheroni constant and Kq is the 
modified Bessel function of the second kind. Therefore, 
c is a weakly dependent function of d for d > a. 

In the case oi W < a/2 additional passages with slopes 
corresponding to noninteger multiples of a/ {2d) open and 



Cl, C2 are modified. Ci is numerically determined, while 
for C2 we have the formula 
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where Nm is the number of those integers in the interval 
[0, m—l] which are relative prime to m. (Note that 
Ni = 1). 
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FIG. 3: Top panel: the return probability P{s) (solid line) and 
its integral F{s) (dashed line) as functions of s/d for W = a/2 
and d — a. In the inset F{s) (solid line) and its asymptotic 
form Fa{s) given by Eq. (O (dashed line) are plotted as func- 
tions of s/d for large s. Bottom panel: P{s) (solid line and 
left y ajcis) for larger s along with the asymptotic expression 
of the return probability Pais) — dFa{s)/ds (dashed line and 
right y axis). 

Using the smooth asymptote of F{s) given by Eq. Q 
and Eqs. Q, an analytic expression can be found for the 
low energy behavior of the DOS: 
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and i?T — M(5/(47r) is the Thouless energy j^. Besides 
the term proportional to the energy (which was also pre- 
dicted in Refs. 0,0), there is an additional logarithmic 
factor. Moreover, in contrast to the above mentioned ref- 
erences the coefficient a depends on W., a and d. However, 
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/3 depends only on W/a. Note that W = a is a special 
case corresponding to the system studied by de Gennes 
and Saint- James Then, C2 — and no logarithmic 
factor appears in the DOS and the result is the same as 
in Ref . iGj 

It is interesting to see the W ^ a limit for which one 
would expect that P(s) has a universal limiting form, 
namely it is only a function of W'^/A. In this case, the 
electron has enough time to explore the whole available 
phase space before escaping (even in an integrable bil- 
liard) and P(s) looses any detailed dependence on the 
geometry of the billiard. Therefore, two billiards with the 
same W'^ /A but different aspect ratio (d/a) should have 
the same path length distribution for small enough W. 
From numerical calculations we found that for — > 0, 
a converges to a value a « 3.62 independent of the as- 
pect ratio, i.e., it becomes a universal constant. Similarly 
we find that when W —^ 0, (3 also tends to a universal 
constant /3 S/tt^. These results are demonstrated in 
Fig. 01 In the inset of the top panel of Fig. 01 the asymp- 
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Top panel: a as 
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function of Wja for dja 



1,2,1/2,1/4,1/8. Bottom panel: /3 as a function of Wja 
obtained from Eq. <8dll using Eqs. 16bl and for C2. The 
inset of the top panel shows the asymptotic DOS na{E) (in 
units of 1/(5) for d = a, W/a = 1,0.5,0.1,0.01 (dot-dashed, 
dotted, long-dashed, and dashed lines, respectively) and the 
result from Ref. 3 (solid line). 

totic DOS na{E) obtained from Eq. Ij8b|l is plotted for 
different values of W along with the result from Ref. 2i 
for the sake of comparison. Our results are in the same 
order of magnitude as that found in Melsen et al. Q and 
Ihra et al. jj]. However, from our analysis it turns out 
that the functional form of the asymptotic DOS is not 
just a linear function but involves a logarithmic factor. 

In conclusion, we have shown that exact quantum me- 
chanical calculations for the integrated DOS of rectan- 
gular Andreev billiards agrees well (for the whole energy 



range below the gap) with that obtained from the Bohr- 
Sommerfeld approximation provided the energy depen- 
dent phase shift is taken into account. From the exact 
analytic form of the asymptote of the integrated return 
probability, we predict a new, logarithmic contribution 
to the DOS on the scale of the Thouless energy. In con- 
trast to earlier results, we show that the DOS at this 
energy range explicitly depends on W, a and d but it has 
a universal limiting form for small enough W. We also 
investigated the case when the superconductor is placed 
anywhere at the side of the rectangle and found that the 
logarithmic contribution in the DOS is generic for rect- 
angular Andreev billiards. 
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